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ON LYMAN-LIMIT SYSTEMS AND THE EVOLUTION OF THE INTERGALACTIC IONIZING BACKGROUND 
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ABSTRACT 

We study the properties of self-shielding intergalactic absorption systems and their implications for 
the ionizing background. We find that cosmological simulations post-processed with detailed radiative 
transfer calculations generally are able to reproduce the observed abundance of Lyman-limit systems, 
and we highlight possible discrepancies between the observations and simulations. This comparison 
tests cosmological simulations at overdensities of ~ 100. Furthermore, we show that the properties 
of Lyman-limit systems in these simulations, in simple semi-analytic arguments, and as suggested 
by recent observations indicate that a small change in the ionizing emissivity of the sources would 
have resulted in a much larger change in the amplitude of the intergalactic H 1-ionizing background 
(with this scaling strengthening with increasing redshift). This strong scaling could explain the rapid 
evolution in the Lya forest transmission observed at z « 6. Our calculations agree with the suggestion 
of simpler models that the comoving ionizing emissivity was constant or even increasing from z = 3 to 
6. Our calculations also provide a more rigorous estimate than in previous studies for the clumping 
factor of intergalactic gas after reionization, which we estimate was «2 — 3atz = 6. 
Subject headings: cosmology: theory — large-scale structure of universe — quasars: absorption lines 



1. INTRODUCTION 

The hydrogen be tween galaxies w as reionized more 
than 12.5 Gyr ago (|Fan et all 120061 ). After this event, 
a largely uniform H i-ionizing radiation background 
pervaded the intergalactic medium (IGM) and kept 
the hydrogen highly ionized everywhere except within 
rare, dense pockets (|Gunn fc Peterson! 119651: iCen et al.l 
19941: iMiralda-Escude et al.lll996t iHernquist et al.l ll99€ 



Haardt fc Madaul ll996UKatz et all 119961 ) . Theampli- 
tude of this background appears to have declined qu ickly 
above a redshift of z — 6 (e.g., iFan et aTJ I2006D and 
to h ave stayed relatively consta nt over 2 < z < 4 
(e.g.. iFaucher-Giguere et aT1l2008[ ). It is of some debate 
whether this decline owed to the overlap stage of cosmo- 
logical reionization or to something more munda ne (e.g., 
IGnedinl[2000tlFan et aLllMl IBecker eTaT]l2007l) . In ad- 
dition, the dense self-shielding systems that remained in 
the IGM, called Lyman-limit systems, provide a win- 
dow into structure formation in a different density regime 
than explored by other large-scale structure probes. 

This paper aims to study intergalactic radiative trans- 
fer and the properties of Lyman-limit systems (systems 
with H I columns of 10 17 2 < N m < 10 19 cm" 2 ). To 
do so, we post-process cosmological simulations (both 
with and without galactic feedback implementations) 
with ionizing radiative transfer. The abundance of 
Lyman-limit systems in the post-processed simulations 
is then compared with measurements of their abun- 
dance from quasar absorption line studies. This com- 
parison tests how cosmological simulations fare at the 
outskirts of galactic halos. Previous attempts to model 
Lyman-limit systems in simulations had reported vari- 
ous levels of agreement with observations and had used 
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much smaller box sizes and particle numbers than are 
presently feasible (iKatz et al.lll996t iGardner et al.ll200H 
IKohler fc Gncdin 2007). In addition, recent observations 
have significantly improved the c onstraints on the abun- 
dance of Lyman-limit system s (|Pro chaska e t al.l 120091 
[20101: ISongaila fc Cowidl20i0l) . 

Our calculations also provide an estimate for the H I- 
ionizing emissivity of star-forming galaxies and quasars 
since, after cosmological reionization, the production rate 
of ionizing photons was in balance with the number 
of absorptions. These absorptions occurred primarily 
within Lyman-limit systems. Our determinations of the 
emissivity using the simulations agree with observations 
at redshifts where it can be measured observationally 
(2 < z < 4), and they also provide a means to estimate 
the emissivity at higher redshifts. We also study the rela- 
tionship between the emissivity of ionizing photons and 
the intensity of the ionizing background. Interestingly, 
we find that small changes in the emissivity at high red- 
shifts could have resulted in much larger changes in the 
ionizing background. This finding has interesting impli- 
cations for the observed trends in this background. 

This paper is organized as follows. Section [5] de- 
scribes the observations and our numerical techniques. 
Section [3] compares our numerical calculations of 
d 2 M/dz cWht - the number of absorption systems per 
unit redshift per H I column (JVhi) - with observa- 
tions. Section 2] studies the relationship between dense 
absorption systems, the ionizing emissivity, and the 
ionizing background. Section [5] considers how physical 
effects that may not be captured in the simulations 
could alter our conclusions. Finally, Appendix [A] 
provides supplementary information regarding the 
properties of the simulated Lyman-limit systems. Our 
calculations assume a flat ACDM cosmology with 
(Cl m ,Cl b ,h,a s ,n s ) = (0.28,0.046,0.70,0.82,0.96) and 
a hydrogen mass fraction of 0.75, consistent with the 
latest Wilkinson Microwave Anisotropy Probe analysis 
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(jKomatsu et al.ll2011l) . 



2. OBSERVATIONS AND MODELS OF H I ABSORBERS 

2.1. Observations 

iProchaska et al.1 ([20091 ) and lProchaska et al.1 ([20101 ) re- 
cently presented measurements of d 2 N ' jdz dNm over the 
interval 10 16 < Nm < 10 19 cm~ 2 . The shaded regions in 
Figure [J (ex cluding the pink r egion) show the constraints 
compiled in Prochask a et al.l ([2010) on 



f(Nm) 



H{z) 



d 2 N 



dzdN m H(0)(l + z) 2 



(1) 



where H(z) is the Hubble parameter. The function 
f(Nm) is independent of redshift if the number den- 
sity and proper size o f the absorbers does not evolve. 
IProchaska et al.l ([2010D inferred that the power-law in- 
dex of /(iVm), — P, must be steeper than —1.7 for 
A/hi < 10 17 5 cm" 2 . This bound constrains the density 
profile of overdensity ~ 100 systems, and we will show 
has interesting implications. In addition, the red error 
bars in the inset i n the left panel o f Figu re [T] show the 
novel constraint of IProc haska ct al. (2009j) on the mean 
free path of 1 Ry photons, A, obtained by stacking spec- 
tra. 

However, the column-density distribution of systems 
wi th 10 16 < Nm < 10 19 cm" 2 is not measured directly 
in IProchaska et all ([20101 ). Rather, it is inferred by in- 
terpolating between the constraints on (1) the number 
density of systems with N m > 10 17 5 cm" 2 , (2) A, (3) 
f{Nm) at columns probed by the Lya forest, and (4) 
measurements at Nm. > 10 19 cm" 2 . The specifics of t his 
interpolation a re described in 



We note that|0'Mcara ct al 



Prochaska etaTI (j2010h . 
(|2007| ) concluded that the 



normalization of the fit resulting in the tan region may 
be be biased high. This study found a large excess of ab- 
sorbers in the smallest Nm. bin, which could owe to spu- 
riously identified systems. When this bin was excluded 
from their analysis, it shifted the best-fit normalization 
of /(AW) down by a factor of ~ 1.5. The pink region in 
Figure [T] assumes a normalization that is reduced by a 
slightly larger factor of 2 (which turns out to be the fac- 
tor needed to agree with our fiducial calculation; Section 
12.21) . For this alternative normalization, the cyan region 
i n Figure U should also be adjusted to enforce continuity. 

iSongaila k, Cowiel ([20101 ) recently measured the abun- 
dance of Lyman-limit systems between < z < 6, which 
extended this measurement to higher redshifts than re- 
ported in previous studies. The magenta points with 
error bars in the inset in Figure [1] show this measure- 
ment, using P = 1.3 as was assumed in their study to 
extrapolate to smaller systems than were measured. Us- 
ing instead /3 = 1.8 would reduce the amplitude of these 
points by a factor of 0.3. 

2.2. Simulations 

Almost all cosmological simulations take the inter- 
galactic gas at z < 6 to be in photoionization equilib- 
rium with a homogeneous radiation field. However, the 
approximation of a uniform ionizing background is not 
appropriate in systems that can self-shield to this radia- 
tion. To allow for self-shielding, we post-process smooth 



particle hydrodynamic s simulations t hat were run with 
the GADGET-2 code ([Springe!! 120051 ) with a specialized 
radiative transfer algorithm. 

Our fiducial simulation uses a 512 3 gas particles and 
the same number of dark matter particles in a box 
of size 20/ft. comoving Mpc (cMpc), resulting in dark 
matter and gas particle masses of 5.4 x 10 6 M© and 
1.1 x 10 6 Mq respectively!! These simulations resolve 
the > 10 9 " 10 Mq halos that are expected to retain their 
gas after reionization with thousands of particles and, 
thus, are the most likely hosts of self-shielded gas. How- 
ever, we have also checked for convergence with regard 
to cosmic variance and resolution by comparing with a 
40//i cMpc simulation run with the same parameters and 
the same number of particles0 

In addition, we also use a 40/ft. cMpc, 2 x 512 3 
simulation with a prescription for galactic winds to 
gauge the importance of such feedback on our re- 
sults. This simulation uses the wind implementation of 
iSpringel fc Hernquistl ([20031 ) with a mass loading factor 
of 2 and with a wind velocity of 680 km s . In this wind 
implementation, twice the mass formed in stars is ejected 
in a biconical wind with approximately the said velocity, 
and the wind is briefly decoupled from the gas such that 
it can escape from the galaxy. This simulation should be 
considered a relatively extreme model for winds based 
on the hea yy suppression of the galaxy mass function it 
yields (see lFaucher-Giguere et al.ll201lD . 

Next, we use the SKID algorithm^ to locate dense re- 
gions in the simulation by linking all particles within 0.5 
mean inter-particle spacings (with the additional crite- 
rion that the baryonic overdensity of associated particles 
is > 2.2) and place each particle group onto a 120 3 rect- 
angular grid for the radiative post-processing. The typi- 
cal size size of this region is a few tenths of a cMpc. We 
have run the same calculations on a 60 3 grid and found 
excellent convergence. We have also verified convergence 
with respect to the linking length choice. 

Each group is illuminated with a power-law inten- 
sity per unit frequency with index —a from 1 to 
4 Ry, and its amplitude is set by the H I photoion- 
ization rate V. In the fiducial calculation, we use 
a = 1 and T_i2 = 0.5, where T_i2 is Y in units 
of 10" 12 s . These choices are in line with mea- 
sure ments at 2 < z < 5 and UV background mod- 
els dBolton fc Haehneltl 120071: lHaardt fc Madaul 119961 : 



iFaucher-Giguere et al.ll2008L l2009h . We find that chang- 
ing a by ±1 alters f(Nm) with respect to the fiducial case 
by < 10% at z = 4, with the largest effect at columns of 
10 18 < Nm < 10 20 cm" 2 . The radiative transfer algo- 
rithm casts six plane-parallel inward fronts of rays along 
the axes of the grid. The rays are attenuated by the 
H I photoelectric optical depth, and the intensity that 

3 Over and above our radiative transfer calculations, the ionizing 
background is turned off on-the-fly in these simulations for regions 
with number densities of hydrogen greater than 10~ 2 cm" 3 in or- 
der to mimic self-shielding as described in Fauchcr-Gigucrc ct al. 
(2010). Wc find that this switch has little impact on our results, 
suppressing /(-/Vjji) by a maximum of 20% at 3 X 10 18 < JVjji < 
3 X 10 20 cm -2 and having no effect at smaller columns. 

4 We find that the normalization of /(A^hi) m the larger box is 
10% higher fairly uniformly for iVjji < 10 20 cm -2 , while it under- 
sho ots by 10% at N m ~ 10 21 cm" 2 . 

5 http : //www-hpec . astro .Washington. edu/tools/skid. html 
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Fig. 1. — Comparison of /(JVhi) in our calculations with observationa l constraints. The shaded r egions show the observ ational constraints 
at a mean redshift of 3.7 fPro chaska et al.H2009l ; lU'Meara et al.H2007l ; IProchaska fc Wolfe 20091; rFrochaska et al.l[201Cf) . Left panel: The 
black thick solid curve is the fiducial calculation, which takes T_i2 = 0.5, the blue dashed curves are the same calculation but with 
r_i2 = 0.25 or 1 (upper and lower curves, respectively), and the black thin solid curve assumes no attenuation of the ionizing background 
(i.e., optically thin radiative transfer). Right panel: The black thick solid curve is the fiducial calculation, the red dotted curve is from the 
simulation with galactic winds, the magenta short-dashed curve uses the alternative recombination prescription, and the blue long-dashed 
curve assumes a gas temperature of 10 4 K in self-shielding regions. The aforementioned curves in both panels are for z = 4, and the 
black thin solid curve in the right panel is the same as the black thick solid but for z = 6. All curves in the right panel are calculated 
for r_i2 = 0.5. The inset in the left panel shows the mean free path of 1 Ry photons i n proper Mpc. The red er ror bars in the inset are 
the [Prochaska ct al. (2009) measurement (2<r errors) and the magenta error bars are the lSongaila jc Cowiel 1 120101) measurement assuming 
/3 = 1.3 (lo" errors). 



reaches each cell sets the local H I fraction under the 
assumption of ionization equilibrium and using the sim- 
ulations' original gas temperatures. The He I fraction is 
assumed to trace the H I fraction with the rest of the 
helium in He II, which is a good approximation in H I 
Lyman- limit systems ([McQuinn fc Switzerl l2010l ). The 
algorithm then iterates to achieve convergence. The as- 
sumptions of a power-law specific intensity and plane- 
parallel sources aligned with the grid allow us to use an- 
alytic expressions for the average photoionization rate in 
a cell as a function of the total H I column to the cell 
and the H I column across the cell itself^ 

3. COMPARISON WITH OBSERVATIONS 

The results of our fiducial calculation at z — 4 are given 
by the black thick solid curve in both panels of Figure [TJ 
This curve has a similar slope at TVhi < 10 175 cm~ 2 to 
what has been inferred observationally (green shaded re- 
gion). The flattening at higher columns, whic h physically 
occurs because of the onset of self-sh ielding (jKatz et al.l 
119961 iZheng fc Miralda-Escu dc 2002), is less prominent 
in the fiducial calculation than in the observations. (The 
effect of this flattening in the fiducial calculation can be 
gleaned by comparing with the black thin solid curve, 
which assumes that the ionizing background is not at- 

6 The ionization front can span less than a proper kpc at relevant 
densities and, thus, is often not resolved in our calculations. We 
find that iterating on the volume-averaged photoionization rate in 
a (uniform density) cell rather than, for example, the cell-centered 
value is critical to converge properly. 
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Fig. 2. — Same curves and shaded regions as in the left panel of 
Figure[T] but the y-axis plots Nm f(Nm) rather than f(Nm). 

tenuated.) However, the fiducial model is consistent at 
2 a with two numbers: the number density of systems 
with Nya > 10 17 ' 5 cm~ 2 at z = 4 and the mean free path 



4 



of ionizing photons (see inset in the left panel of Fig. [T]), 
the two constraints that principally determine the green 
and cyan regions. For a finer comparison, see Figure [2j 
which instead plots Nm /(-/Vhi) for the same quantities 
as in the left panel of Figure [1] 

The level of agreement with observation s is an im- 
provem ent over previous numerical studies. iKatz et all 
(|1996() and lGardner et al.l (|2001f l found that Lyman-limit 
systems were underproduced by an order of magnitude 
in their cosmological simulations, likely because of the 
lo wer resolution o f thei r calculations. The simulations 
m IGardner et all pool were run m a 11 h^ 1 cMpc box 
with 2 x 64 3 particles and had particle masses that were 
roughly two orders of magnitude larger than in our sim- 
ul ations. Our results are in better agreement with those 
of iKohler fc Gnedinl (|2007l ). who, for a suitable choice of 
F, was able to reproduce the abundance of Lyman-limit 
systems in their 4 h~ Y cMpc and 8 cMpc simulations 
with 128 3 dark matter particles and the same number 
of Lagrangian gas elements. However, our /(iVm) has 
more structure than th e power law-like form found in 
IKohler fe Gnedinl (|2007t ). Our calculations qualitatively 
agree with those of Altay et al. (2011, submitted con- 
currently with this study), which used simulations with 
similar box sizes and particle numbers to those discussed 
here. 

The form of f(Nm) in our calculations depends on 
several factors. The thin solid curve in the right panel 
in Figure [T] is f(Ngj) at z = 6 for the fiducial back- 
ground model. This curve is similar to the z = 4 
curve, but is steeper at Nm < 10 18 cm -2 . In ad- 
dition if the background is varied so that T_ 12 = 
0.25 or 1 (which spans the range of most estimates; 
iFaucher-Giguere et al.l |2008|) . /(Ahi) becomes respec- 
tively the upper and lower blue dashed curves in the left 
panel of Figure [TJ The fiducial calculation does not in- 
corporate winds that recycle gas from within galaxies. 
The red dotted curve in the right panel is /(-/Vhi) calcu- 
lated from the same ionizing background model as the 
fiducial model curve, but with a simulation that incor- 
porates galactic winds as described in Section |2"T2"1 The 
abundance of Lyman-limit systems is increased slightly 
in the case with winds compared to our fiducial calcula- 
tion because winds result in more gas surrounding halos. 
We have also looked at a simulation that assumes half of 
the wind velocity and half the mass loading factor as the 
aforementioned case and found surprisingly similar re- 
sults compared to this model with more extreme winds. 

Our calculations assume Case A recombinations. This 
neglects the reabsorption of ionizing photons produced 
by recombinations to the ground state. The magenta 
short-dashed curve in the right panel in Figure [T] lo- 
cally accounts for this reabsorption by using as + (a a — 
as)r'/r as the recombination coefficient in a grid cell, 
where V is the local (attenuated) value of the photoion- 
ization rate, and a a and o.b are the Case A and B re- 
combination coefficients. This interpolates between Case 
A and Case B recombinations in optically thin and thick 
regions, respectively. This improved prescription sup- 
presses the abundance of Ahi ~ 10 19 cm~ 2 systems at 
the 10% level. A countering effect to neglecting reabsorp- 
tion is that the simulations' temperatures are likely over- 
estimated in self-shielding regions owing to suppressed 



atomic cooling f q5.1[) . The blue long-dashed curve is 
for a toy temperature model where the temperature is 
set to f = [10 4 + (T - 10 4 ) r'/r] K, where T is the 
original temperature in the simulation. Over relevant 
columns, To ~ 15 — 20 x 10 4 K in the simulations so that 
T < To in this model. The primary difference between 
this model and the fiducial model is that /(iVm) is in- 
creased at Nm ~ 10 19 5 cm -2 . While this model likely 
overcompensates for the effect of this additional cooling, 
it does not account for the impact of lower temperatures 
on the gas dynamics ( £15.11) . 

4. THE CASE FOR A STRONG RELATIONSHIP BETWEEN 
EMISSIVITY AND T 

We now focus on the implications of the functional 
form of /(-/Vhi). In particular, we show that its form 
implies that small changes in the emissivity of the sources 
can result in much larger changes in T. 



4.1. The \Miralda-Escude et~ali (MM) Model 

We first describe the semi-analytic model of 
iMiralda-Escude et al.l [20001 (henceforth MHR), which is 
useful for interpreting the observations and numerical 
calculations. The MHR model assumes that there is 
a critical density Aj above which systems self-shield 
from the ionizing background and remain neutral. We 
will show that our calculations support this assump- 
tion. Given Aj, the volume filling fraction of neutral 
gas, Fy, and the global recombination rate density, 
71, can be calculated from just the volume-weighted 
gas density probability distribution -P(A), a quantity 
easily measured from cosmological simulations. For 

example, for isothermal gas 1Z oc J Q dAA 2 P(A). We 
will subsequently use that 1Z was in balance with the 
ionizing emissivity of the sources, e, which is a good 
approximation at z > 2 and after reionization (e.g., 
IMadau et aflll999l) . 

This model further assumes that the mean free path 
at 1 Ry is given by A = ao Fy(Ai)~ 2 / 3 , where ao is an 
adjustable parameter. This relation for A would hold if 
the comoving number density and shape of self-shielding 
regions were independent of A^ . The primary conclusion 
that we derive from the MHR model would be strength- 
ened if the number density decreased or the shape be- 
came more spherical with Aj, as one would anticipate. 

Intuition is gained by approximating P(A) as a power 
law with index —7 over the densities that contri bute most 
of the recombinations (|Furlanetto fc Ohll2005D . In this 
approximation, the MHR model yields the relations 



A= a 



-1 -2/3 



dAP(A) 



A, 



A 



2( 7 -l)/3 



e= K = a A n e n H / dAA 2 P(A)ocA 



(2) 
(3) 



r _ <tl l (3^-3 + Q!) c ^ a (7- 7 )/3 
3 + a 



^ OC f V(2-« 



(4) 

where nx = nxA is the number density in species X, 
<tll is the photoionization cross section at 1 Ry, — j8 is 
the power law index of f(Nm) at iVni ~ 10 17 cm -2 , and 
we have assumed that the gas is isothermal (which holds 
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Fig. 3.— A 3 times the PDF of the gas density, P(A), at the 
specified redshift. The curves in the panels on the left are A 3 -P(A) 
from published gas density PDFs and from our calculations. The 
curves labeled "RT" in the panels on the right only include the 
contribution from photoionized gas as described in the text, which 
results in a cutoff at high densities. The key in the upper panels 
applies to both the upper and lower panels. 

approximately in our simulations at relevant densities). 
Equation ^ is valid for a > 0, and it assumes that 
cosmological expansion is unimportant (valid at z > 2) 
and that the mean free path of an ion izing photon with 
energ y E equals A \Ej\ Ryl 3/3 - 3 (e.g.. lHaardt fc Madaul 
119961 iMiralda-Escudel 12003) . The final proportionality 
in equation (j4|) assumes spherically symmetric power- 
law density profiles such that 7 can be related to the 
(3 that applies for optically thin systems via the rela- 
tions f(N m )dN m oc lixrdr, P(A) dA cx 4irr 2 dr, and 
riHi oc A (i.e., there is an annulus at radius r with col- 
umn ./Vhi and a shell at r with density A). In the previous 
section, we found that /3 « 1.8 for optically thin systems 
in both the simulations and in observations. Equation 
((I]) indicates that such a (3 results in a strong relation- 
ship between e and T. 

The power-law of the density PDF, 7, is more funda- 
mental for determining this relationship than the power- 
law of f(Nni), f3. At high densities, the gas density PDF 
in the simulations can be approximated with the power- 
law 2.5 < 7 < 3, which also implies a strong scaling 
between e and T (eqn. |4}. The left panels in Figure 
[3] show A 3 P(A) at z = 4 (top right panel) and z — 6 
(bottom right panel) The dashed red curve is the fit to 
P( A) in MHR. The b lue dotted curve is the fit to more 
recent simulations bv iBolton fc Beckerl [20091 (BB). The 

7 The small break at A = 10 in the curves in Figure \3\ originates 
because we are patching together P(A) from a calculation that uses 
a linking length of 1 mean interparticle spacing with our fiducial 
calculation that uses a linking length of 0.5 in order to cover the 
plotted range in A. 
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Fig. 4. — Top Panel: Number of H 1-ionizing photons emitted 
per hydrogen atom per Gyr as a function of redshift. The blue 
circles are the fiducial calculation, and the red squares are the 
calculation with winds. The different points correspond to r_i2 = 
0.1, 0.25, 0.5, and 1, with the point size increasing with T_i2. The 
error bars bracket the observational constraints on r_i2, assuming 
the fiducial model mapping from r_i2 to e. Bottom panel: The 
power-law index that relates e to T from the fiducial (blue solid 
curve) and wind (red dotted curve) simulations. These curves are 
computed from the ratio of the emissivity between the T_i2 = 0.25 
and 0.5 calculations. 

black solid curve is A 3 P(A) in the simulation without 
winds, and the magenta dot-dashed curve is the simu- 
lation with strong winds. The simulation with winds 
has more gas around galaxies, which increases A 3 P(A) 
somewhat relative to the fiducial simulation^ 

Interestingly, all of the P(A) yield 7 > 2.5, where —7 
is the power-law index of P(A). For the MHR fit to 
P(A), 7 « 2.5 and equation ((J]) implies that T oc e" 
where n = 30 Thus, the intergalactic H 1 photoioniza- 
tion rate scales as the cube of the emissivity. However, 
n is extremely sensitive to 7 in this equation. The P(A) 
in our simulations and in those of BB show 7 = 2.5 — 3 
at relevant A (with 7 increasing with redshift), which 
formally imply n = 3 — 00 given equation (H]). Therefore, 
the flatness of A 3 P(A) indicates a strong relationship 
between e and T. 

How justified is the MHR assumption of a critical den- 
sity at which gas becomes neutral? The curves labeled 
"RT" in the right panels of Figure [3] are the PDF of ion- 
ized gas in our calculations, or more specifically A 3 P(A) 
multiplied by Tnni/[a.A n^]. The latter factor is equal 
to the square of the ionized fraction in the absence of 
collisional ionizations. The recombination rate is pro- 

8 The different P(A) curves should not be expected to agree 
for a few reasons. These include that they were computed from 
simulations which had different thermal histories and because the 
MHR fit was to a simulation using an outdated cosmology and with 
a rigid paramctrization. 

9 For power-law density profiles, 7 = 2.5 corresponds to isother- 
mal spheres with A oc r -2 , and 7 = 3 corresponds to A oc r . 
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portional to the logarithmic integral over the RT curves. 
The RT curves follow the density PDF curves before an 
abrupt cutoff, which supports the MHR model assump- 
tion of a characteristic A;. 

4.2. Numerical Calculations 

We can measure the power-law relationship between 
e and T directly in the simulations by calculating the 
recombination rate of the gas (which is in balance 
with e) as a function of rPn Recombinations that re- 
sult from collisional ionizations are omitted by count- 
ing photo-ionizations and using that this is equal to the 
photoionization-induced recombination rate. The indi- 
vidual points in Figure |4] are calculated directly from the 
fiducial simulation (blue circles) and the simulation with 
winds (red squares). Each set of points at fixed z corre- 
sponds to this calculation for T_i 2 = 0.1, 0.25, 0.5, and 
1, with the point size increasing with T_i2E3 

Our numerical calculations indicate that T should scale 
strongly with emissivity. At z = 3, we find n = 2.2 — 2.4 
in the simulations (bottom panel, Fig. |4|. At z — 4, n = 
2.5 — 2.8, which is slightly weaker than the MHR model 
prediction for 7 = 2.5 of n — 3. Section [5T2l provides an 
explanation for this weaker dependence. At z = 6, we 
find an even stronger scaling, with n = 3.8 — 4.5. Thus, a 
small change in e, particularly at z — 6, could result in a 
much larger change in T. For example, if e changes from 
2.1 to 2.5 photons per hydrogen atom per Gyr at z — 6, 
then r_i2 evolves from 0.1 to 0.25. We argue in Section 
[6] that this behavior could expl ain the quick e volution 
in r that is observed at z w 6 (Fa n et alJ 12006ft . which 
many studies have interpreted as indicating the overlap 
stage of reionization. 

These calculations also constrain the intergalactic 
emissivity history. The error bars in Figure [4] bracket 
the blue po ints' T_ 12 values that are consistent with mea- 
surem ents (jFaucher-Giguere et al.1 [20081 : ICalverlev et al.1 
2011). In particular, they bracket r_i 2 = {0.6 — 
1.4, 0.25-1, 0.2-0.7, 0.07-0.3} at z = {3,4,5,6}. Each 
error bar generously spans a factor of w 4 in r_i2, which 
is much larger than the error quoted in many studies and 
spans the range of most previous estimates. However, 
since e depends weakly on T in these calculations, the 
uncertainty in the inferred emissivity is much smaller. 

If r_i2 evolves from s» 1 at z = 3 to » 0.1 at z = 6 as 
observations suggest, this evolution results in a slight in- 
crease in the comoving emissivity towards higher redshift 
in our calculations. The primary difference between our 
weakly increasing comoving emissivity history and the 
weakly decreasing history of Bolton & Hachnelt (2007) 
[calculated using constraints on T and the A calculated 
with the MHR model] is that our mean free path val- 
ues decrease more rapidly with increasing redshift. At 
z — 6, a comoving emissivity of approximately two ion- 
izing photons per hydrogen atom in the age of the Uni- 

10 Or, we could equivalently calculate the rate of absorptions 
given f(Nm), and we have verified that this yields identical results. 

11 Our group-finding method is 50% complete at A = 5 and 
more complete at higher A. We assume the MHR fit for P(A) 
and isothermal 20, 000 K gas to add the contribution to TZ from 
lower densities. This correction is largest in the calculation at 
2 = 6 and r_i2 = 0.1, where this gas contributes 25% of the 
total recombinations. Taking lower temperatures for this gas would 
result in more recombinations and would increase n somewhat. 



verse (1 Gyr) is needed in our calculations to be consis- 
tent with the constrain ts on r_i2- This results i n the 
same dilemma posed in iBolton &; Haehnehl (|2007f ) that 
this comoving emissivity was only sufficient to reionizc 
the Universe by z = 6 if e were constant or increasing 
to higher redshift. The increasing nature of our inferred 
comoving emissivity history with increasing redshift at 
z > 3 may point to the solution. 

These calculations also provide an estimate for the 
intergalactic clumping factor (the enhancement in the 
recombination rate over a homogeneous universe) since 
this number is equal to e/\aAn 2 H \. Studies that esti- 
mate the ionizing emissivity from observations of z > 6 
galaxies need to know this uncertain factor to quantify 
whet her this emissivity can keep the Universe ionized 
(e.g. JMadau et al.|[i~99a iRobertson et al.ll2010D . Unfor- 
tunately, there is a history of numerical simulations using 
very large subgrid clumping factors of ~ 30 to compen- 
sate for recombinations in "unresolved" structures. (The 
number 30 originated from counting gas that was ionized 
by loc al sources in the simulations of lGnedin fc Ostrikerl 
(|1997f ) in the clumping factor tabulation. The dumpi- 
ness of this gas, and hence the higher rate of recombina- 
tions in it, is better accounted for as a n escape fraction 
of ionizing photons; IKohler et al J 120071 ) Recent studies 
have calculated the clumping factor from simulations a s 
a function of A;, finding 1 - 10 (e.g. JPawlik et al.N2009h . 
Because our calculations determine the ionization struc- 
ture of self-shielding regions, we can estimate the clump- 
ing factor directly rather than leave it as a function of A, . 
In fact, the logarithmic integral over the RT curves in the 
right panel in Figure [3] equals the clumping factor. We 
find a clumping factor enhancement of 2.4 — 2.9 at z = 6 
for the observationally allowed range of T_i2 = 0.1 — 0.25 
and using the simulation without winds (but see tj5.4[) . 
We find similar values using the simulation with winds. 
Here we define the clumping factor as the enhancement 
in the recombination rate over a homogenous universe 
with a gas temperature of 10 4 K. 

5. ADDITIONAL CONSIDERATIONS 

The results of this study rely on the density profile 
of A ~ Aj systems being adequately simulated and the 
post-processing approximation for the radiative transfer. 
Here we address briefly how physical processes that are 
neglected in our calculations could alter the results. 

5.1. Cooling Physics 

The cooling of self-shielding gas is not properly cap- 
tured in cosmological simulations without self-consistent 
ioniz ing radiative transfer (e.g., iFaucher-Giguere et al.1 
2010). In particular, the collisional cooling rate is pro- 
portional to the H 1 fraction times the electron density, 
which is generally underestimated in the simulations at 
Njjx > 10 17 cm~ 2 since the temperature of the gas is 
computed assuming a homogeneous ionizing background. 
(This error is more important than the cosmological sim- 
ulations' underestimate of the photohcating rate in self- 
shielding, photo-ionized gasQ) Because collisional cool- 
ing is very efficient for gas with temperature greater than 

12 It is only when systems are not highly photoionized - TVht > 
10 19 ' 5 cm -2 in our calculations - that the photoheating rate is 
overestimated in the simulations. 
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re 20, 000 K and inefficient at lower temperatures, under- 
estimating the amount of collisional cooling may not lead 
to a serious error. We find that if we set the tempera- 
ture to 10 4 K in self-shielding regions in post-processing 
as described in Section [31 this operation flattens /(iVm) 
somewhat but does not affect n. Lower temperatures and 
the correct ionization levels would also affect the distri- 
bution of the gas through their effect on the gas pressure. 
See footnote 3 for a comparison that is relevant for esti- 
mating the significance of these effects. 

In ad dition, our simu l ations do not include metal line 
cooling. iProchter et al.l (|2010| ) measured the metallicity 
of a Lyman-limit system and a few super Lyman-limit 
systems at z re 3.5 to be Z re 10~ 17 Zq, similar to 
the average metallicity of damped Lya systems. This 
metallicity is insufficient for metals to be the dominant 
coolant at T > 20,00 K, which requires Z > lO" 1 Z Q 
(jWiersma et al.l[2009| ). However, this metallicity is suffi- 
cient to allow ionized gas to cool to somewhat lower tem- 
peratures than in the simulations, where (T) re 20, 000 K 
at relevant columns. 

5.2. Collisional Ionizations 

Feedback from galaxies could result in the gas within 
Lyman-limit systems being shock heated, driving it into 
a collisional rather than photoionization equilibrium. If 
the Lyman-limit systems were in collisional ionization 
equilibrium, their cross section would be insensitive to F, 
so that r would scale linearly with the emissivity. In our 
simulations with and without galactic winds, collisional 
ionizations are unimportant for A ~ 100 gas at z = 6, 
but become more important with decreasing redshift and, 
in fact, set the ionization state for half of A ~ 100 gas 
at z = 3. In fact, we find that collisional ionizations play 
some role in making our calculations' n smaller at z < 4 
than the predictions based on equation (|4]). 

As a final note, if /? in reality were a s steep at TVht 5, 
10 175 cm" 2 as the observations of iProchaska et al.l 
(2010) suggest, collisional ionization is likely the only 
means to significantly weaken the dependence of F on e 
from the scaling in equation 

5.3. Proximity Effects 

While the extragalactic ionizing background almost 
certainly dominates in systems with JVhi < 10 172 cm -2 
(|Miralda-Escudeir2005D . local sources of radiatio n could 
be im portant for higher column-density systems (|Schavd 
2006). Our calculations do not account for their contri- 
bution. The analytic model in §4.1l can be generalized to 
allow for a fraction of the gas at each density, /(A), to 
be ionized by local sources. In this model, the T depen- 
dence of the emissivity would only be weakened relative 
to our previous findings if -P(A) [1 — /(A)] decreased less 
quickly with A than P(A). However, the opposite seems 
more likely to hold, i.e., that local sources preferentially 
ionize the denser regions near the sites of star formation. 

5.4. The Reionization History 

Reionization is only crudely included in our simula- 
tions, with the intergalactic gas being reionized homoge- 
neously at z re 10 in these calculat ions. However, reion- 
ization could end as late as z k 5.5 (|Mesingerl l2010). Our 
results may be affected if reionization occurred at lower 



redshifts than is assumed in the simulations. The relax- 
ation timescale after a system was reionized is roughly 
equal to the Jeans length divided by the sound speed, 
or H(z)^ 1 A -1 / 2 . It takes Az re A" 1 / 2 (1 + z) after a 
region was reionized for the gas to have relaxed. For 
the A ~ 100 systems of interest, this corresponds to 
Az ~ 0.1 (1 + z). After a large-scale region was reionized, 
the local P(A) would decrease, starting at the highest 
densities and evolving to lower densities. This inside-out 
process would likely steepen the scaling between T and 

— 1/2 

e a time interval of Az > A i (1 + z) after a region 
was reionized. These considerations would also result in 
our calculations underestimating the dumpiness of the 
ionized gas just after reionization, and thereby the emis- 
sivity, at fixed r. 

6. DISCUSSION 

The abundance of Lyman-limit systems in our fidu- 
cial calculation is consistent with observations at the 2 a 
level. We find that adding feedback from galactic winds 
or lowering T improves the agreement. This comparison 
tests whether our cosmological simulations reproduce the 
properties of gas at the outskirts of galactic halos. How- 
ever, our simulations may underproduce the number of 
Nm ~ 10 19 cm -2 systems compared with observations. 
If real, this discrepancy could indicate that the gas is 
cooler and, as a result, denser in these systems than in 
the simulations owing to additional cooling, or that the 
impact of feedback is more prominent than in the con- 
sidered wind model. 

Both observations and numerical simulations suggest 
a steep slope for f(Nm) with (3 re 1.8 at Ahi < 
lO 17 5 cm~ 2 . We showed that this steepness implies a 
strong relationship between the extragalactic H i-ionizing 
background and the ionizing emissivity of the sources. 
At z = 3, our calculations find that T scales as e to 
the power of 2.2 — 2.4. The near constancy of T that 
is measured at z w 2 - 4 (|Bolton fc Haehnelti 120071: 
IFaucher-Giguere et al.ll2008fl puts stringent requirements 
on the evolution of the total comoving ionizing emis- 
sivity owing to this strong scaling, requiring it to be 
roughly constant. This is an interesting finding since 
the ionizing background is believed to transition from 
being dominated by quasars to being dominated by 
star-forming galaxies over this redshift interval (e.g., 
IFaucher-Giguere et al1l2008H . 

The scaling strengthens with redshift in our calcula- 
tions, with r proportional to e to the 3.8 — 4.5 power 
at z = 6. A strong scaling is related to the claim that 
the IGM clumpin g factor after reioni zation was small, in- 
dependent of Ai (jPawlik et al.ll2009h . since both require 
7 k 3 for the high-density power- law index of the gas 
density distribution. Interestingly, iFan et ail (|2006t ) de- 
tected a sharp decrease in the Lya forest tr ansmission 
at z re 6. If interpreted as a change in T, IFan et aT] 
(2006) estimated that a factor of at least 2 change was 
required over Az re 0.5. It does not seem plausible that 
the emissivity of galaxies changed by such a large factor 
over this interval, which constitutes 7% of the Hubble 
time at z = 6. Our simulations find that a more mod- 
est re 20% change in the comoving emissivity at z re 6 
over Az re 0.5 would result in a factor of 2 change in T. 
Estimates for the variance in T(z) between widely sep- 



arated regions suggest that this effect could operate on 
this timescale (Mcsinge r fe Furla nctto 20091). 

The quick evolution in T from the mechanism discussed 
here is not necessarily associated with the end of reion- 
ization, defined as the epoch when neutral patches were 
present in the low-density IGM so that A,- ~ 1. In all 
of our calculations, A; 3> 1, with Aj ~ 50 at z f» 6 and 
for r = 10~ 13 s _1 - approximately what is measured at 
this z using the Lya and Ly/3 forests. In addition, the 
mean free path of 1 Ry photons is 8 cMpc for this case. 
This is still larger than the 0.5 (1.1) cMpc average sep- 
arations of galaxies residing in > 10 8 (> 10 9 ) M Q halos 
in the assumed cosmology. Therefore, it does not satisfy 
the simplest criterion for reionization that the mean free 
path was sho rter than the separation of the anticipated 
sources (e.g., iMiralda-Escude et al.ll2000l )FI 

The leading theory for explaining the rapid evolution 
in r had been that it owed to the overlap of large-scale 
cosmological H n regions (e.g.. lGnedinll2000T ). However, 
it appears unlikely that overlap can produce such rapid 
evolution in T for two reasons. First, models predict large 
temporal variance between when different ~ 30 cMpc lo- 
cations in the IGM were ionized, with t he majority of the 
gas ionized over an interval of Az > 3 (|Furlanetto et al.l 
12004 iMcQuinn etaH I2007D . Second, by the time of 
overlap, it is likely that T was regulated by absorptions 
in Lyman-limit systems rath er than in diffuse H I gas 
(|Furlanetto fc Mesingerll2009D . If this were the case, the 
overlap process itself could not mediate a sharp jump in 
A (and hence T), and fast temporal evolution in T could 
only owe to the type of effect discussed hereJ3 

Our fiducial calculations favor the comoving emissiv- 
ity of H i-ionizing photons to have increased slightly 
with redshift from z — 3 to 6. This trend is consis- 
tent with the conclusion that the emissivity had to be no 
more than weakly decreasing wit h redshift at z > 3 to 
be able to reionize the Universe ([Miralda-Escudel [20031 : 



iBolton fc Haehnelti [20071 : iFaucher-Giguere et al.l [20081) . 

The improvement in this work over these previous stud- 
ies is a self-consistent model for X(T), which is necessary 
to infer the emissivity at z > 4. We also showed that 
this conclusion is robust to the considerable uncertainty 
in r at z « 6. A growing emissivity with redshift could 
arise from the efficiency of escape of ionizing photons in- 
creasing with redshift or as a result of thermal feedback 
on low-mass galaxies. 

In contrast with this r e sult, si mulations of reion- 
ization (e g.. Illiev et al.l [2001 iTrac fc Cenl 120071 : 
iZahn etaH 20071 ) generally prescribe an emissivity 
history that is strongly decr easing with increasing 
redshift (jChoudhurv et al.l feOOQ). In addition, previous 
radiative transfer simulations of reionization in large 
enough box sizes to study the statistics of this process 
(> 10 comoving Mpc; e.g. . Illiev et alJl2006t ITrac fc Cer] 
120071: IZahn et al.ll2007T ) did not resolve the Lyman-limit 
systems and, thus, could not capture the effects dis- 
cussed in this paper. To capture these systems at z > 6 
necessitates resolving the ionization structure of gas 
parcels with A ~ 100. While this would require a large 
increase in resolution beyond that achieved in these 
previous simulations of reionization, a first step would 
be to incorporate Lyman-limit systems with a subgrid 
model based on the calculations described here. 
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Fig. 5. — Average properties of absorbers as a function of iVjii in the fiducial calculations. The solid blue curves are for 2 = 4, and the 
red dotted curves are for 2 = 6. Each curve is calculated by weighting by jihi along a skewer and then averaging skewers at fixed Nm. The 
top panel shows the average of T_i2 within absorbers. For comparison, the two black thin curves in the top panel are the same quantities 
but for the two analytic slab models described in Appendix [A] The middle panel shows the average H I fraction. The bottom panel shows 
the average density in units of the cosmic mean, and the one-sided error bars shown for the z = 4 case are th e standard deviation in A. 
For comparison, the black thin dashed lines in the bottom panel represent the analytic model of Schayc (2001) at z = 6 and 2 = 4 (lower 
and upper lines, respectively). 
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APPENDIX 

A. PROPERTIES OF LYMAN-LIMIT SYSTEMS 

The panels in Figure [5] show the average of the H I photoionization rate, the fraction of hydrogen in H I, and the 
density. Each property is tabulated by weighing by rim along a skewer through a radiatively post-processed group 
and then averaging over all skewers through all groups at fixed -/Vhi- In each panel, the thick blue curve is the fiducial 
model at z = 4, and the red dotted curve is the same but at z = 6. 

The top panel shows the average of T_i2 as a function of TVhi- For comparison, the black thin dashed and dotted 
curves in this panel are respectively the average of T_i2 in a slab with column JVhi exposed to a planar source radiating 
orthogonal to the slab's plane and with (1) a monochromatic ionizing background at 1 Ry (thin dashed curve) and 
(2) the same ionizing background as in the fiducial model (thin dotted curve). Both curves assume T_i2 = 0.25 to 
be incident on the slab in order to separate these curves from the simulation curves which took T_i2 = 0.5. The 
attenuation of T in slab model (2) is almost identical to that in the simulations. Both decrease as iV^j 1 at large 
columns. 

The middle panel shows the average H I fraction as a function of Nm- Systems with Nm < 10 19 5 cm -2 are highly 
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ionized in our calculatio ns. This is consist ent with observations, w hich find that Lyman-limit systems are highly 
ionized (Prochaska 1991) and the results of iKohler fc Gne din (2003). It occurs because (1) harder photons require 
larger iVm to be attenuated (a 4 Ry photon has o ptical depth unity for iVm = 1 x 10 19 cm" 2 ) and (2) even T_i2 <C 1 
can be sufficient to keep the gas highly ionized (|Zheng fc Miralda-Escudl 120021 ). In our calculations and those of 
IKohler fc Gnedinl ((20071 ). systems with 10 17 < N m < 10 19 cm 2 are significantly more ionized than those in the study 
of lAltav et al.l (|2010t ). For example, systems with AThi = 10 18 cm -2 in lAltav et al.l (|2010D had a median H I fraction 
of 0.1 — 0.3 versus 10~ 2 here. 

The bottom panel in Figure [5] shows the average of A. The density of absorbers increases as « before plateauing 
over 10 18 < A^hi < 10 20 cm~ 2 owing to self-shielding. This plateau occurs at A ~ 500 for z — 4 and at A ~ 100 for 
z — 6. The one-sided error bars represent the standard deviation in the density at z = 4. Prior to self-shielding, 
there is a tight relationship between density and A/hi because the abso rbers' size is roughly the Jeans length and 
their ionization state is set by the amplitude of the ionizing background jSchaye 2001). Once self-shielding becomes 
significant, this relationship is broken. As a result, the scatter becomes much larger at A^hi 10 17 cm~ 2 . This scatter 
decreases again when the systems become neutral at A^hi > 10 20 cm -2 as a similar relationship is again established. 
Lastly, if we had plotted the average proper density rather than the density in units of the mean in the bottom panel, 
the z = 4 and z = 6 curves would almost overlap. 

The black thin dashed lines in the bottom panel of Figure [5] are the predicted relationship in the model of iSchavel 
(|2001f ) between A and A^i for optically thin systems with T_ 12 = 0.5 and 2 x 10 4 K gas (although, the temperature 
dependence in this model is weak). This model just assumes photoionization equilibrium with the average T and that 

2/3 

the absorbers width is the Jeans length, predicting the scaling A oc A^ HI . We find that this simple model underprcdicts 
A by a factor of « 2 at A^hi ~ 10 16 cm" 2 , with this difference decreasing with increasing Nm. However, this model 
should not be expected to capture the normalization at better than the factor of 2 level, and, more importantly, we 
find that it is successful at capturing the trend with A^hi- 



